% Created 2018-11-08 Thu 17:29
% Intended LaTeX compiler: xelatex
\documentclass[11pt]{article}
\usepackage{graphicx}
\usepackage{grffile}
\usepackage{longtable}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage[normalem]{ulem}
\usepackage{amsmath}
\usepackage{textcomp}
\usepackage{amssymb}
\usepackage{capt-of}
\usepackage{hyperref}
\usepackage{ctex}
\author{\href{mailto:pengwenbin7@126.com}{pengwenbin7}}
\date{\today}
\title{机器学习概述与线性回归}
\hypersetup{
 pdfauthor={\href{mailto:pengwenbin7@126.com}{pengwenbin7}},
 pdftitle={机器学习概述与线性回归},
 pdfkeywords={},
 pdfsubject={},
 pdfcreator={Emacs 26.1 (Org mode 9.1.9)}, 
 pdflang={English}}
\begin{document}

\maketitle
\tableofcontents


\section{概述}
\label{sec:org3b62958}
机器学习三要素: 模型、策略、算法

实现机器学习的步骤：
\begin{enumerate}
\item 得到一个有限的训练数据集合
\item 确定包含所有可能的模型假设空间，即学习模型的集合
\item 确定模型选择的准则，即学习的策略
\item 实现求解最优模型的算法，即学习的算法
\item 通过学习方法选择最优的模型
\item 利用学习的最优模型对数据进行预测或分析
\end{enumerate}

常用数据预处理方法：
\begin{itemize}
\item 标准化, \(\mu\) 为平均值， \(\sigma\) 为标准差:\(\overrightarrow{x} \gets (\overrightarrow{x} - \mu) / \sigma\)
\item 缩放， \(\lambda\) 为常数: \(\overrightarrow{x} \gets \lambda \overrightarrow{x}\)
\item 对数变换:\(\overrightarrow{x} \gets log(\overrightarrow{x})\)
\end{itemize}

常用损失函数：
\begin{itemize}
\item 0-1损失函数: \$L(Y, f(X)) = \begin{cases}
1, \& Y \(\neq\) f(X) \\
0, \& Y = f(X) \nonumber
\end{itemize}
\end{cases}\$
\begin{itemize}
\item 平方损失函数: \(L(Y, f(X)) = (Y - f(X)) ^ 2\)
\item 绝对损失函数: \(L(Y, f(X)) = |Y - f(X)|\)
\item 对数损失函数: \(L(Y, f(X)) = -logP(Y | X)\)
\end{itemize}

损失函数的期望称为 \textbf{风险函数(risk function)} 或 \textbf{期望损失(expected loss)}
$$
R_{exp} (f) = E_{p}[L(Y, f(X))] = \int_{x, y}L(y, f(x))P(x, y)dxdy
$$
离散情况下, \textbf{经验风险(empirical risk)}
$$
R_{emp} (f) = \frac{1}{N} \sum_{i = 1} ^{N} L(y_i , f(x_i))
$$
根据大数定律，当样本容量 N 趋于无穷时， R\(_{\text{emp}}\) 趋向于 R\(_{\text{exp}}\) 。按照 \textbf{经验风险最小化} 求解最优模型:
$$
argmin \frac{1}{N} \sum_{i = 1} ^{N} L(y_i , f(x_i))
$$
为了防止 \textbf{过拟合(over-fitting)}, 加上正则化项，结构风险最小化
$$
R_{srm} (f) = \frac{1}{N} \sum_{i = 1} ^{N} L(y_i , f(x_i)) + \lambda J(f)
$$

\section{线性回归}
\label{sec:org83aa70d}
\href{https://en.wikipedia.org/wiki/Linear\_regression}{wiki:线性回归}

\subsection{一元回归}
\label{sec:org960b19f}
\url{https://upload.wikimedia.org/wikipedia/commons/thumb/3/3a/Linear\_regression.svg/438px-Linear\_regression.svg.png}

对于上图，蓝点 \((x_i, y_i)\) 是已知的数据，求红线 \(y = \alpha x + \beta\)，
称 x 为 \emph{特征} ， 用于预测当 x 为 \(x_0\) 时 y 的 \emph{估计值} \(\hat{y}_0\)，
设 \emph{训练集} （用于拟合的数据集合）的样本数量为 m ，定义 \emph{代价函数} (损失函数)： 
$$
C(\alpha, \beta) = \frac{1}{2m}\sum_{i = 1}^{m}{(\hat{y}_{i} - y_{i})^2}
= \frac{1}{2m}\sum_{i = 1}^{m}{(\alpha x_{i} + \beta - y_{i})^2}
$$

\begin{quote}
代价 = 均方差 / 2, 除 2 的意义在于消去导数的系数2，方便计算。
\end{quote}

当代价函数取得较小的值时，这时的 (\(\alpha\), \(\beta\)) 是较好拟合。
接下来可以使用 \emph{梯度下降} 、 \emph{最小二乘法} 等方法求解。
显然，当只有1组数据 (x, y) 时，即为平面上一个点。过这一点，可以做任意条直线，
代价都是0。对于多元回归也是类似的， 即训练数据的数目 m 不小于特征值的数目 n。

\subsection{多元回归}
\label{sec:org9cd9eb8}
相对于一元回归，多元回归就是 x 变成 \(\overrightarrow{x}\) 向量。以经典的\href{http://archive.ics.uci.edu/ml/machine-learning-databases/housing/}{波士顿房价} 为例。
\begin{center}
\label{tab:org2a3a7db}
\begin{tabular}{rrrlr}
0.00632 & 18.00 & 2.310 & \dots{} & 24.00\\
0.02731 & 0.00 & 7.070 & \dots{} & 21.60\\
0.02729 & 0.00 & 7.070 & \dots{} & 34.70\\
\dots{} & \dots{} & \dots{} & \dots{} & \dots{}\\
0.84054 & 0.00 & 8.140 & \dots{} & 13.90\\
\end{tabular}
\end{center}

上表中，每一行中前面的列表示输入，最后一列表示输出，
列的意义为城镇犯罪率、房间数、黑人比例等特征，这里不去考虑它们的实际意义。
如果仅考虑第一列和最后一列的关系，即为 \emph{一元回归} ，
如果考虑前面列和最后一列的关系，即为 \emph{多元回归} 。

\begin{quote}
记特征数为 n (列数 - 1)，行数为 m 。
以 x\(^{\text{(i)}}_{\text{j}}\) 表示输入数据第 i 行第 j 项。 y\(^{\text{(i)}}\) 表示第 i 个输出数据。
\end{quote}

对于第 i 组数据，设拟合函数为 \(h_{\theta}(x^{(i)})\), 估计值为 \(\hat{y}^{(i)}\)，有:
$$
\hat{y}^{(i)} = h_{\theta}(x^{(i)}) = \theta_0 + \theta_1 x^{(i)}_1 + \theta_2 x^{(i)}_2 + \cdots + \theta_n x^{(i)}_n \tag{1}
$$

为了计算方便，设 x\(^{\text{(i)}}_{\text{0}}\) = 1 ，称为 \emph{偏置单元} ，则 (1) 式表示为:
$$
\hat{y}^{(i)} = \sum_{j = 0}^{n}\theta_j x^{(i)}_j
$$

第 i 组的代价为: \(\frac{(\hat{y}^{(i)} - y^{(i)})^2}{2}\)。
所以，总体代价为:

\begin{align*}
C(\theta) &= \frac{1}{2m} \sum_{i = 1}^{m}(\hat{y}^{(i)} - y^{(i)})^2  \\
&= \frac{1}{2m} \sum_{i = 1}^{m} [(\sum_{j = 0}^{n}\theta_j x^{(i)}_j) - y^{(i)}]^2  \\
&= \frac{1}{2m}
\begin{pmatrix}
(\theta_0 x^{(1)}_0  + \theta_1 x^{(1)}_1 + \theta_2 x^{(1)}_2 + \cdots + \theta_n x^{(1)}_n - y^{(1)})^2 + {} \\
(\theta_0 x^{(2)}_0  + \theta_1 x^{(2)}_1 + \theta_2 x^{(2)}_2 + \cdots + \theta_n x^{(2)}_n - y^{(2)})^2 + {} \\
\vdots \\
(\theta_0 x^{(m)}_0  + \theta_1 x^{(m)}_1 + \theta_2 x^{(m)}_2 + \cdots + \theta_n x^{(m)}_n - y^{(m)})^2
\end{pmatrix} \tag{2}
\end{align*}

求 C(\(\theta\)) 关于于 \(\theta_{\text{j}}\) 偏导

\begin{align*}
\frac{\partial C}{\partial \theta_{j}} &= \frac{1}{m} \sum_{i = 1}^{m}(\hat{y}^{(i)} - y^{(i)}) \frac{\partial \hat{y}^{(i)}}{\partial \theta_{j}}
\nonumber \\
 &= \frac{1}{m} \sum_{i = 1}^{m}(\hat{y}^{(i)} - y^{(i)}) x^{(i)}_j
\nonumber \\
&= \frac{1}{m} \sum_{i = 1}^{m}(h_{\theta}(x^{(i)}) - y^{(i)}) x^{(i)}_j \tag{3}
\end{align*}

最后用\href{https://en.wikipedia.org/wiki/Gradient\_descent}{梯度下降法}求解 \(\theta_{\text{j}}\) \\

$$
\text{iterate}: \\
\theta_j \gets \theta_j - \frac{\alpha}{m} \frac{\partial C}{\partial \theta_{j}}
$$

其中 \(\alpha\) 称为 \emph{学习率} ,通常取一个小随机数，且每一次迭代使用的 \(\alpha\) 可以不一样。

考虑结构风险最小化， 代价函数

\begin{align*}
J(\theta) &= C(\theta) + \frac{\lambda}{2} \sum_{j = 1}^{n + 1} \theta_{j} \\
&= \frac{1}{2m} \sum_{i = 1}^m (\hat{y}^{(i)} - y^{(i)}) ^ 2 + \frac{\lambda}{2} \sum_{j = 1}^{n + 1} \theta_{j}^2 \tag{4}
\end{align*}

(4) 中 \(\frac{\lambda}{2} \sum_{j = 1}^{n + 1} \theta_{j}^2\)
称为正规化项， \(\lambda\) 称为 \textbf{正规化参数} ， 写成 \(\frac{\lambda}{2m}\),
下标 j 从 0 或 1 开始，都对运算结果无影响。由 (4) 式也可以看出，当 \(\lambda\) 过大时，
代价函数对 \(\theta_{j} (j > 0)\) 的惩罚过大，此时拟合函数接近于 \(y = \theta_0\)
表示的直线，形成 \textbf{欠拟合} 。
对应的梯度为

$$
\frac{\partial J}{\partial \theta_j} = 
\frac{1}{m} \sum_{i = 1}^{m}(h_{\theta}(x^{(i)} - y^{(i)}) x_{j}^{(i)} + \lambda \theta_j \tag{5}
$$

\subsection{向量化}
\label{sec:orge749205}
首先，将 \(x\), \(y\), \(\theta\) 写成矩阵的形式

\begin{align*}

\overrightarrow{X} &=
\begin{bmatrix}
x^{(1)}_{0} & x^{(1)}_{1} & \cdots & x^{(1)}_{n} \\
x^{(2)}_{0} & x^{(2)}_{1} & \cdots & x^{(2)}_{n} \\
\vdots      & \vdots       & \ddots & \vdots      \\
x^{(m)}_{0} & x^{(m)}_{1} & \cdots & x^{(m)}_{n} \\
\end{bmatrix} \\

\overrightarrow{y} &= 
\begin{bmatrix}
y^{(1)} \\
y^{(2)} \\
\vdots   \\
y^{(m)}
\end{bmatrix} \\

\overrightarrow{\theta} &=
\begin{bmatrix}
\theta_0 \\
\theta_1 \\
\theta_2 \\
\vdots \\
\theta_n
\end{bmatrix}

\end{align*}

\subsubsection{梯度下降}
\label{sec:orgee0c60f}
(1), (2), (3) 式对应的向量表达如下：

\begin{equation*}

\overrightarrow{\hat{y}} = h_{\overrightarrow{\theta}}(\overrightarrow{X}) = \overrightarrow{X} \overrightarrow{\theta} \\

C(\overrightarrow{\theta}) = \frac{1}{2m} (\overrightarrow{\hat{y}} - \overrightarrow{y})^{T} (\overrightarrow{\hat{y}} - \overrightarrow{y}) \\ 

\frac{\partial C}{\partial \overrightarrow{\theta}} = \frac{1}{m} \overrightarrow{X}^{T} (\overrightarrow{X} \overrightarrow{\theta} - \overrightarrow{y})  \\

\end{equation*}

梯度下降求 \(\theta\)
$$
\text{iterate}: \\
\overrightarrow{\theta} \gets \overrightarrow{\theta} - \frac{\alpha}{m} \frac{\partial C}{\partial \overrightarrow{\theta}}
$$

若考虑结构风险最小化，令

\begin{align*}

\overrightarrow{L} = \lambda \begin{bmatrix}
0 & 0 & 0 &\cdots & 0 \\
0 & 1 & 0 & \cdots & 0 \\
0 & 0 & 1 & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & 1 \\
\end{bmatrix}

\end{align*}

\begin{verbatim}
L = lambda * eye(length(theta));
L(1) = 0;
\end{verbatim}

有

\begin{align*}

J(\overrightarrow{\theta}) = \frac{1}{2m} (\overrightarrow{\hat{y}} 
- \overrightarrow{y}) ^T (\overrightarrow{\hat{y}} - \overrightarrow{y}) 
+ \frac{1}{2} {\overrightarrow{\theta}^T \overrightarrow{L} \overrightarrow{\theta}} \\

\frac{\partial J}{\partial \overrightarrow{\theta}} = 
\frac{1}{m} \overrightarrow{X}^{T} (\overrightarrow{X} \overrightarrow{\theta} - \overrightarrow{y})
+ \overrightarrow{L} \overrightarrow{\theta}

\end{align*}

\subsubsection{最小二乘法}
\label{sec:org72f76b8}

由 \(\overrightarrow{y} = \overrightarrow{X} \overrightarrow{\theta}\) 可以直接算出
\(\overrightarrow{\theta} = (\overrightarrow{X}^{T} \overrightarrow{X}) ^{-1} \overrightarrow{X}^{T} \overrightarrow{y}\)

\begin{quote}


\(\overrightarrow{X}\) 为 \(m \times (n + 1)\) 矩阵，仅当 \(m = (n + 1)\) 时，其逆矩阵存在，而 \(\overrightarrow{X}^{T} \overrightarrow{X}\) 是 \(n + 1\) 阶
方阵，其逆存在，由逆矩阵的性质 \(\overrightarrow{A} \overrightarrow{A}^{-1} = \overrightarrow{E}\) 及矩阵乘法结合律 \((\overrightarrow{A} \overrightarrow{B}) \overrightarrow{C} = \overrightarrow{A} (\overrightarrow{B} \overrightarrow{C})\)，
对 \(\overrightarrow{y} = \overrightarrow{X} \overrightarrow{\theta}\) 两边同时右乘 \((\overrightarrow{X}^{T} \overrightarrow{X})^{-1} \overrightarrow{X}^{T}\) ，
得到
\begin{align*}
(\overrightarrow{X}^{T} \overrightarrow{X})^{-1} \overrightarrow{X}^{T} \overrightarrow{y} &= (\overrightarrow{X}^{T} \overrightarrow{X}) ^ {-1} \overrightarrow{X}^{T} \overrightarrow{X} \overrightarrow{\theta} \\
&= (\overrightarrow{X}^{T} \overrightarrow{X})^{-1} (\overrightarrow{X}^{T} \overrightarrow{X}) \overrightarrow{\theta} \\
&= \overrightarrow{E} \overrightarrow{\theta} \\
&= \overrightarrow{\theta}
\end{align*}
\end{quote}


\section{代码实现}
\label{sec:org1073e71}
实现载入数据、一元回归样本、最小二乘法和梯度下降的工具函数 \texttt{tools.py}
\begin{verbatim}
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import os
import urllib
import numpy as np


def loadData():
    """ 载入数据 """
    # 下载数据文件
    dataFile = 'housing.data'
    if not os.path.isfile(dataFile):
	urllib.request.urlretrieve(
	    'http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data',
	    'housing.data')

    # 使用 numpy.loadtxt 载入数据
    data = np.loadtxt(dataFile)
    # 观察到数据中存在从某一位置开始全为0的列，做平滑处理
    data = data + 1e-9
    return data


def oneSamples(m=50):
    """ 预处理用于一元回归的样本数据 """
    data = loadData()
    # 取前m组数据
    x = data[0:m, 0]
    y = data[0:m, -1]

    # 排序
    index = np.argsort(x)
    x = x[index]
    y = y[index]

    # 标准化: (原始值 - 平均值) / 标准差
    x = (x - np.mean(x)) / np.std(x)
    y = (y - np.mean(y)) / np.std(y)

    # 数据转矩阵，不习惯二维数组
    X = np.matrix(x).transpose()
    Y = np.matrix(y).transpose()

    # 加入偏置单元 X <- [ 1, X ]
    X = np.hstack([np.ones((m, 1)), X])

    return (X, Y)


def leastSquare(X, Y):
    """ 最小二乘法计算参数矩阵 theta
    Parameters
    ----------
	X: 输入矩阵, m x (n + 1)
	Y: 输入矩阵, m x 1
    Returns
    -------
	theta: 参数矩阵, (n + 1) x 1
    """
    return np.linalg.pinv(X.transpose() * X) * X.transpose() * Y


def gradientDescent(theta,
		    X,
		    y,
		    alpha=1e-3,
		    lparameter=0,
		    maxIter=1000,
		    stopCost=1e-6):
    """ 
    梯度下降实现
    Parameters
    ----------
	X: 包含偏执单元的输入样本矩阵，大小为 m x (n + 1), m 为样本个数， n 为特征个数
	y: 输出样本矩阵, 大小为 m x 1
	alpha: 学习率，梯度下降的速率
	lparameter: 正规化参数
	maxIter: 最大迭代次数
	stopCost: 当损失小于此值时，提前终止迭代并返回结果
    Returns
    -------
	theta: 待求参数, (n + 1) x 1 的矩阵
	J: 损失
	counter: 迭代次数
    """
    counter = 0
    m = X.shape[0]
    J = np.zeros(maxIter)
    L = lparameter * np.eye(theta.shape[0])
    L[0][0] = 0

    while counter < maxIter:
	h = X * theta - y
	cost = 1 / (2 * m) * h.transpose() * h + \
	    theta.transpose() * L * theta / 2
	J[counter] = cost
	counter += 1
	if cost < stopCost:
	    break
	else:
	    theta = theta - alpha / m * X.transpose() * (
		X * theta - y) + L * theta

    J = J[0:counter]

    return (theta, J, counter)
\end{verbatim}

观察函数图像，调节参数 \texttt{ex.py}
\begin{verbatim}
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import tools

# 取前50组数据方便看图
(X, Y) = tools.oneSamples(50)

# 显示样本数据
plt.plot(X[:, 1], Y, color='k', label='sample')
plt.legend()

# 最小二乘法求解 theta
theta = tools.leastSquare(X, Y)

# 求估计值
hatY = X * theta

# 画出图像
plt.plot(X[:, 1], hatY, color='r', label='least squares')
plt.legend()

# 使用梯度下降求解
th = np.matrix([[0], [0]])

# 传入不同的 alpha, 观察函数图像
(theta, J, counter) = gradientDescent(th, X, Y, alpha=1e-3)
plt.plot(X[:, 1], X * theta, color='g', label='grad-0.001')
plt.legend()

(theta, J, counter) = gradientDescent(th, X, Y, alpha=0.01)
plt.plot(X[:, 1], X * theta, label='grad-0.01')
plt.legend()

(theta, J, counter) = gradientDescent(th, X, Y, alpha=0.1)
plt.plot(X[:, 1], X * theta, label='grad-0.1')
plt.legend()

# 保存图片，非必须
plt.savefig('plot.png')
\end{verbatim}
运行以上代码，得到如下图片:
\begin{center}
\includegraphics[width=.9\linewidth]{plot.png}
\end{center}

通过观察，当 \(\alpha\) 为 0.1 时，拟合结果和最小二乘法结果基本重合。
当 \(\alpha\) 小于 0.01 时，因为迭代次数不够，尚未很好的拟合。
另外，可以通过画出不同学习率下，损失函数的图像，寻找合适的学习率 \texttt{costReg.py}
\begin{verbatim}
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import tools

(X, Y) = tools.oneSamples()
th = np.matrix([[0], [0]])

(theta, J, ci) = tools.gradientDescent(th, X, Y, alpha=2)
plt.plot(J)
plt.savefig('a-2-J.png')
plt.close()

(theta, J, ci) = tools.gradientDescent(th, X, Y, alpha=0.1)
plt.plot(J, label='a=0.1')
plt.legend()

(theta, J, ci) = tools.gradientDescent(th, X, Y, alpha=0.01)
plt.plot(J, label='a=0.01')
plt.legend()

(theta, J, ci) = tools.gradientDescent(th, X, Y, alpha=0.001)
plt.plot(J, label='a=0.001')
plt.legend()
plt.savefig('a-0.1-0.0001-J.png')
\end{verbatim}

当 \(\alpha\) = 2 时，损失函数如下：
\begin{center}
\includegraphics[width=.9\linewidth]{a-2-J.png}
\end{center}

可以看到，函数不能收敛，说明学习率过大。

当 \(\alpha\) 为 0.001， 0.01， 0.1 时
\begin{center}
\includegraphics[width=.9\linewidth]{a-0.1-0.0001-J.png}
\end{center}

函数均收敛，选取收敛最快的 0.1 作为学习率。
\begin{quote}
因为特征值经过标准化处理后，待求的参数本身和 [0, 0] 非常接近，所以这里可以非常快地收敛。

数组和矩阵的转化是不必须的，看着习惯一点。
\end{quote}

\section{调用API}
\label{sec:org19ac791}
\subsection{使用 sklearn 随机梯度下降拟合}
\label{sec:org9db5a86}
\texttt{use-api.py}
\begin{verbatim}
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# 参考: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDRegressor.html

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import SGDRegressor
import tools

clf = SGDRegressor(max_iter=1000)
(X, Y) = tools.oneSamples()
Y = Y.A1

clf.fit(X, Y)
hatY = clf.predict(X)

plt.plot(X[:, 1], Y, label='sample')
plt.legend()
plt.plot(X[:, 1], hatY, label='sgd')
plt.legend()
plt.savefig('sgd.png')
\end{verbatim}
得到拟合图像如下:
\begin{center}
\includegraphics[width=.9\linewidth]{sgd.png}
\end{center}

\begin{quote}
matlab/octave 的 fit 函数在多项式拟合时尤其好用。
\end{quote}

\section{泛化能力}
\label{sec:org42417be}
波士顿房价原始数据有506组，使用交叉验证的方式，前400组作为训练数据，后106组作为测试数据。
比较 \texttt{tool.leastSquare}, \texttt{tools.gradientDescent} 和 \texttt{sklearn.linear\_model.SGDRegressor} 的预测结果。
\begin{verbatim}
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import tools
from sklearn.linear_model import SGDRegressor

# 载入数据
data = tools.loadData()

# 按列对数据做标准化处理
oriD = np.matrix(data)
(M, N) = oriD.shape
D = np.matrix(np.zeros((M, N)))
for i in range(0, N):
    mu = np.mean(oriD[:, i])
    sigma = np.std(oriD[:, i])
    D[:, i] = (oriD[:, i] - mu) / sigma

# 取前 400 组数据作为训练数据
m = 400
stdX = D[0:m, 0:-1]
Y = D[0:m, -1:]
(m, n) = stdX.shape

# 加入偏执单元
X = np.hstack([np.ones((m, 1)), stdX])

# 最小二乘法求解
thetaLS = tools.leastSquare(X, Y)

# 初始化 theta
th = np.matrix(np.zeros((n + 1, 1)))

# 梯度下降求解
(thetaGD, J, ci) = tools.gradientDescent(th, X, Y, alpha=0.1, maxIter=1000)

# 训练 sgd
clf = SGDRegressor(max_iter=1000)
clf.fit(stdX, Y.A1)

# 取 400 以后的数据作为测试数据
stdX = D[m:, 0:-1]
Y = D[m:, -1:]
(m, n) = stdX.shape

# 加入偏执单元
X = np.hstack([np.ones((m, 1)), stdX])

# 最小二乘
hatLS = X * thetaLS
# 梯度下降
hatGD = X * thetaGD
# sgd
hatSGD = clf.predict(stdX)

plt.plot(Y, label='ori', marker='p')
plt.legend()

plt.plot(hatLS, label='ls', marker='o')
plt.legend()

plt.plot(hatGD, label='gd', marker='+')
plt.legend()

plt.plot(hatSGD, label='sgd', marker='x')
plt.legend()

plt.savefig('generalize.png')
\end{verbatim}
在交互式环境下得到以下结果:
\begin{center}
\includegraphics[width=.9\linewidth]{ge-sc.png}
\end{center}

观察到3条拟合曲线基本重合，可以通过左下角的放大镜按钮局部放大，看到实际是三条曲线。
说明函数的实现是可行的，但和测试数据存在一定偏差，可以对不同列使用不同的预处理方式，
使用高阶拟合、加入正规化项等方式以达到更好的拟合。

\section{更多}
\label{sec:orgace6691}
\subsection{向量化的意义}
\label{sec:orga527e45}
按矩阵乘法的定义计算n阶方阵乘法

$$
\overrightarrow{B} = \overrightarrow{A} \overrightarrow{A} \\
b_{i, j} = \sum_k^{n} a_{ik} a_{kj}
$$

时间复杂度是 \(O(n^{3})\)。

测试环境为 i3-3225, 8G。以下代码计算 1000 阶方阵的乘积，比较使用向量化和不使用向量化的差别。
\begin{verbatim}
clear;close;clc;

n = 1000;
A = ones(n, n);
B = zeros(n, n);

display('iterator:');
tic
for k = 1: n
    for l = 1: n
	for m = 1: n
	B(k, l) = B(k, l) + A(k, m) * A(l, m);
	end
    end
end
toc

display('vector:');
tic
C = A * A;
toc
\end{verbatim}

输出结果为:
\begin{verse}
iterator:\\
Elapsed time is 43.067627 seconds.\\
vector:\\
Elapsed time is 0.044712 seconds.\\
\end{verse}
可以看出，运行时间相差4个数量级。
对于机器学习问题，不实现向量化，往往是没有意义的。

\subsection{改进}
\label{sec:org9976e5e}
\end{document}
